Polyfem Examples¶
This is a jupyter notebook. The "real" notebook can be found here.
Some imports for plotting
import plotly.offline as plotly
import plotly.graph_objs as go
import plotly.figure_factory as ff
#Necessary for the notebook
plotly.init_notebook_mode(connected=True)
algebra
import numpy as np
stuff for the animation
import ipywidgets as widgets
from ipywidgets import interact
and finallypolyfempy
import polyfempy as pf
Plotting utility¶
This code is not particularly interesting.
It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.
def create_plot_mesh_and_layout(vertices, connectivity, function):
x = vertices[:,0]
y = vertices[:,1]
if vertices.shape[1] == 3:
z = vertices[:,2]
else:
z = np.zeros(x.shape, dtype=x.dtype)
if connectivity.shape[1] == 3:
f = connectivity
else:
#Convert a tet-mesh into face based triangles.
f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
for i in range(len(connectivity)):
f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])
mesh = go.Mesh3d(x=x, y=y, z=z,
i=f[:,0], j=f[:,1], k=f[:,2],
intensity=function, flatshading=function is not None,
colorscale='Viridis',
contour=dict(show=True, color='#fff'),
hoverinfo="all")
layout = go.Layout(scene=go.layout.Scene(
aspectmode='data',
xaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
),
yaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
),
zaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
)
))
return mesh, layout
def plot(vertices, connectivity, function, camera=None):
mesh, layout = create_plot_mesh_and_layout(vertices, connectivity, function)
if camera is not None:
layout.scene.camera = camera
fig = go.Figure(data=[mesh], layout=layout)
plotly.iplot(fig)
Creates a quad mesh of n_pts x n_pts in the form of a regular grid
def create_quad_mesh(n_pts):
extend = np.linspace(0,1,n_pts)
x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
pts = np.column_stack((x.ravel(), y.ravel()))
faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)
index = 0
for i in range(n_pts-1):
for j in range(n_pts-1):
faces[index, :] = np.array([
j + i * n_pts,
j+1 + i * n_pts,
j+1 + (i+1) * n_pts,
j + (i+1) * n_pts
])
index = index + 1
return pts, faces
Set the mesh path
mesh_path = "plane_hole.obj"
create settings
settings = pf.Settings()
pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$
settings.discr_order = 1
normalize the mesh to be in $[0,1]^2$
settings.normalize_mesh = True
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
Next we setup the problem
problem = pf.GenericTensor()
sideset 1 has zero displacement in $x$
problem.set_displacement(1, [0, 0], [True, False])
sideset 4 has zero displacement in $y$
problem.set_displacement(4, [0, 0], [False, True])
sideset 3 has a force of [100, 0] applied
problem.set_force(3, [100, 0])
fianally set the problem
settings.set_problem(problem)
Solve!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh_from_path(mesh_path)
solver.solve()
Get the solution
[pts, tets, disp] = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
top_camera = dict(
up=dict(x=0, y=1, z=0),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=0, z=0.8)
)
plot(vertices, tets, mises, top_camera)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call
mises = solver.get_sampled_mises()
plot(vertices, tets, mises, top_camera)